############################################################
#' Spatially explicit individually based model for prey
#' grazing on grass.
#'
#' Author: Jim Graham
#' Date: 3/29/2018
############################################################


library('RColorBrewer')

############################################################
# Prey Functions
############################################################
############################################################
#' Function to create the prey initially
#' @param MaxPrey Maximum number of prey in the model
#' @param NumPrey Number of prey to start the model with
############################################################
PreyInitialize=function(MaxPrey,NumPrey)
{
  xs=rep(0,MaxPrey)
  ys=rep(0,MaxPrey)
  health=rep(0,MaxPrey)
  repotime=rep(0,MaxPrey)
  
  # randomly place the prey on the grass
  
  Index=1
  while (Index<=NumPrey)
  {
    xs[Index]=runif(1, min = 0, max = 1)
    ys[Index]=runif(1, min = 0, max = 1)
    health[Index]=runif(1, min = 0, max = 1) # start health at 1
    repotime[Index]=runif(1, min = 0, max = TimeToRepo) #time since reproduction
    
    Index=Index+1
  }
  ThePrey=data.frame(xs,ys,health,repotime)
  
  return(ThePrey)
}
############################################################
#' Function to move the prey on each cycle
#' @param MaxPrey Maximum number of prey in the model
############################################################
PreyMove=function(MaxPrey)
{
  Index=1
  while (Index<=MaxPrey)
  {
    if (ThePrey$health[Index]>0)
    {
      # randomly move the prey
      X=ThePrey$xs[Index]+rnorm(1,sd=MoveSD)
      #print(ThePrey$xs)
      Y=ThePrey$ys[Index]+rnorm(1,sd=MoveSD)
      
      if (X<0) X=0
      if (X>1) X=1
      if (Y<0) Y=0
      if (Y>1) Y=1
      
      # make sure they stay on the field
      ThePrey$xs[Index]<<-X
      ThePrey$ys[Index]<<-Y
    }
    
    Index=Index+1
  }
}
############################################################
#' Function to have the prey feed on the grass for each cycle
#' @param MaxPrey Maximum number of prey in the model
############################################################
PreyFeed=function(MaxPrey)
{
  NumRows=nrow(TheGrass)
  NumCols=ncol(TheGrass)
  
  # have each prey forage for food
  Index=1
  while (Index<=MaxPrey)
  {
    if (ThePrey$health[Index]>0)
    {
      # make sure we have an x and y coordinate that is valid
      x=as.integer(ThePrey$xs[Index]*NumCols)+1
      if (x>NumCols) x=NumCols
      if (x<1) x=1
      #print(x)
      y=as.integer(ThePrey$ys[Index]*NumRows)+1
      if (y>NumRows) y=NumRows
      if (y<1) y=1
      #print(y)
      
      # get the available food
      AvailableFood=TheGrass[x,y]
      #print(AvailableFood)
      
      # update what happend when they tried to forage
      if (AvailableFood>0) # there is food available
      {
        AmountEaten=MouthFul
        
        AvailableFood=AvailableFood-AmountEaten
        
        if (AvailableFood<=0) # prey ate all the available food
        {
          AmountEaten=-AvailableFood
          AvailableFood=0
        }
        # add the amount eaten to health
        #print(AmountEaten)
        Health=ThePrey$health[Index]+AmountEaten
        if (Health>1) Health=1
        ThePrey$health[Index]<<-Health
        
        #print(AvailableFood)
        TheGrass[x,y]<<-AvailableFood
      }
    }
    Index=Index+1
  }
}
#######################################################
#' Function to determine which prey die when their health
#' goes to zero
#' @param MaxPrey Maximum number of prey in the model
#######################################################
PreyDie=function(MaxPrey)
{
  Index=1
  while (Index<=MaxPrey)
  {
    if (ThePrey$health[Index]>0)
    {
      ThePrey$health[Index]<<-ThePrey$health[Index]-HealthDedecution
      
      if (ThePrey$health[Index]<=0)
      {
        ThePrey$xs[Index]<<-0
        ThePrey$ys[Index]<<-0
        ThePrey$health[Index]<<-0
        ThePrey$repotime[Index]<<-0
        
        #message(cat("Death, NumPrey1=",NumPrey))
        NumPrey<<-NumPrey-1
        #message(cat("Death, NumPrey2=",NumPrey))
      }
      
    }
    
    Index=Index+1
  }
}
#######################################################
#' Function to have prey repoduce
#' @param MaxPrey Maximum number of prey in the model
#######################################################
PreyRepo=function(MaxPrey)
{
  # check for agents that are ready to reproduce
  Index=1
  while (Index<=MaxPrey)
  {
    if (ThePrey$health[Index]>0)
    {
      if (ThePrey$repotime[Index]>=TimeToRepo) # time to give birth
      {
        ThePrey$repotime[Index]<<-0
        
        BirthIndex=1
        while ((BirthIndex<=MaxPrey)&&(ThePrey$health[BirthIndex]!=0))
        {
          BirthIndex=BirthIndex+1
        }
        if (BirthIndex<=MaxPrey)
        {
          # recreate all the vectors in the data frame with new ones!
          # this is something R is not good at
          ThePrey$ys[BirthIndex]<<-runif(1, min = 0, max = 1)
          ThePrey$xs[BirthIndex]<<-runif(1, min = 0, max = 1)
          ThePrey$health[BirthIndex]<<-1
          ThePrey$repotime[BirthIndex]<<-0
          
          NumPrey<<-NumPrey+1
          #message(sprintf("Birth, BirthIndex=%d, NumPrey=%d",BirthIndex,NumPrey))
        }
      }
      else
      {
        # move all the repotimes ahead by one unit
        ThePrey$repotime[Index]<<-ThePrey$repotime[Index]+1
        
      }
    }
    Index=Index+1
  }
}
#######################################################
# Grass functions
#######################################################
#######################################################
#' Function to have the grass grow with each cycle
#' @param MaxPrey Maximum number of prey in the model
#######################################################
GrassGrows=function()
{
  NumRows=nrow(TheGrass)
  NumCols=ncol(TheGrass)
  
  X=1
  while (X<=NumCols)
  {
    Y=1
    while (Y<=NumRows)
    {
      if (TheGrass[X,Y]<1) TheGrass[X,Y]<<-TheGrass[X,Y]+GrassIncrease
      
      Y=Y+1
    }
    X=X+1
  }
  
}
#######################################################
# Utility functions
#######################################################
#######################################################
#' Function to update the plots for the current data.
#' @param TheGrass The matrix with values for the grass
#' @param TheColors A color ramp for the grass
#' @param ThePrey The data table with data for the prey
#######################################################
UpdatePlots=function(TheGrass,TheColors,ThePrey)
{
  # update the plots
  par(mfcol=c(1, 2) )
  
  image(TheGrass,col=rf)
  points(ThePrey,pch = 21:25,bg=TheColors[1],cex=1) # pch is type, cex is size magnifier
  
  plot(Population)
  
  # have R sleep so we can see the graphs
  Sys.sleep(0.2)
}
############################################################
# Global Variables
############################################################

# grass parameters
GrassIncrease=0.001 # Amount the grass increases on each cycle
NumCols=20
NumRows=20

# prey parameters
MaxPrey=1000 # maximum number of prey we can have at one time
NumPrey=20 # number of prey created at the start
MouthFul=0.2 # amount of grass prey eat on each cycle
MoveSD=0.01 # stdev of movement on each cycle (0 to 1)
HealthDedecution=0.1 # amount health declines on each cycle
TimeToRepo=200 # number of cycle suntil prey reproduce

# number cycles to complete
NumIterations=10000

############################################################
# Main Code
############################################################

#Setup a palette of shades of green for the grass
rf <- palette(brewer.pal(n = 9, name = "Greens"))  # make colors
TheColors=heat.colors(10)

# setup the grass for the prey to feed on
z=runif(NumCols*NumRows, min = 0, max = 1)
TheGrass=matrix(z,ncol=NumCols,nrow=NumRows)

#initilaize the prey array
ThePrey=PreyInitialize(MaxPrey,NumPrey)

# Initialize the population vector to track population over time
Population=c(NumPrey)

# repeat the simulation for the specified number of interations
a=0
while (a<NumIterations)
{
  if (a%%10==0)
  {
    UpdatePlots(TheGrass,TheColors,ThePrey)
  }
  # have the prey feed on the grass
  PreyFeed(MaxPrey)
  
  # have the prey move around
  PreyMove(MaxPrey)
  
  # kill off the prey that are weak
  PreyDie(MaxPrey)
  
  # kill off the prey that are weak
  PreyRepo(MaxPrey)
  
  # Have the grass grow back
  GrassGrows()
  
  # update the statistics
  Population=c(Population,NumPrey)
  a=a+1
}